Los modelos jerárquicos involucran varios parámetros de tal manera que las creencias de unos de los parámetros dependen de manera significativa de los valores de otros parámetros de manera jerárquica.

Para ilustrar esta dependencia consideraremos un ejemplo. Consideremos el caso en el que tenemos varias monedas acuñadas en la misma casa de monedas, es razonable pensar que una fábrica sesgada a águilas tenderá a producir monedas con sesgo hacia águilas. La estimación del sesgo de una moneda depende de la estimación del sesgo de la fábrica que a su vez está influido por los datos de todas las monedas. Veremos que la estructura de dependenica a lo largo de los parámetros generan estimaciones mejor informadas para todos los parámetros.

Si pensamos únicamente en dos monedas que provienen de la misma casa de moneda tenemos:

  1. Creencias iniciales de los posibles valores de los parámetros (sesgos de las monedas).

  2. Creencias iniciales de la dependencia de los parámetros por provenir de la misma fábrica.

Cuando observamos lanzamientos de las monedas actualizamos nuestras creencias relativas a los sesgos de las monedas y también actualizamos nuestras creencias acerca de la dependencia de los parámetros.

Ejemplo: Modelo jerárquico de una moneda

Recordemos el caso de lanzamientos de una moneda. Le asignamos una distribución inicial beta, recordemos que la distribución beta tienen dos parámetors \(a\) y \(b\):

\[p(\theta)=\frac{1}{B(a,b)}\theta^{a-1}(1-\theta)^{b-1}\]

Con el fin de hacer los parámetros más intuitivos los podemos expresar en términos de la media \(\mu\) de la distribución inicial y el tamaño de muestra \(K\), el cual refleja la variabilidad de la distribución, esto es, la confianza que tenemos en la media. Los parámetros correspondientes en la distribución beta son:

\[a=\mu K, \quad b = (1-\mu)K\]

Ahora introducimos el modelo jerárquico. En lugar de especificar un valor particular para \(\mu\), consideramos que \(\mu\) puede tomar distintos valores (entre 0 y 1), y definimos una distribución de probabilidad sobre esos valores. Podemos pensar que esta distribución describe nuestra incertidumbre acerca de la construcción de la máquina que manufactura las monedas.

Veremos posteriormente que en el caso de más de una moneda el modelo permite que cada moneda tenga un sesgo distinto, pero ambas tenderán a tener un sesgo cercano a \(\mu\), algunas aleatoriamente tendrán un valor de \(\theta\) mayor a \(\mu\) y otras menor. Entre más grande \(K\) mayor será la consistencia de la acuñadora y los valores \(\theta\) serán más cercanos a \(\mu\). Si observamos varios lanzamientos de una moneda tendremos información tanto de \(\theta\) como de \(\mu\).

Para hacer un análisis bayesiano aún nos hace falta definir la distribución inicial sobre el parámetro \(\mu\); utilicemos una distribución Beta:

\[p(\mu)=Beta(\mu|A_{\mu}, B_{\mu})\]

donde \(A_{\mu}\) y \(B_{\mu}\) se conocen como hiperparámetros y son constantes. En este caso estamos consideramos que \(\mu\) se ubica típicamente cerca de la media \(A_{\mu}/(A_{\mu} + B_{\mu})\) y tomaremos \(K\) como constante.

Recordemos también que en el ejemplo de una moneda teníamos que la función de verosimilitud es una distribución Bernoulli para los datos:

\[p(x|\theta) = \theta^x(1-\theta)^{1-x}\]

Si utilizamos las iniciales Beta para \(\mu\) y \(\theta\) como discutimos arriba, solo nos resta aplicar la regla de Bayes con nuestros dos parámetros desconocidos \(\mu\) y \(\theta\):

\[p(\theta, \mu|x)=\frac{p(x|\theta,\mu)p(\theta,\mu)}{p(x)}\]

Hay dos aspectos a considerar en el problema:

  1. La verosimilitud no depende de \(\mu\) por lo que \[p(x|\theta, \mu)=p(x|\theta)\]

  2. La distribución inicial en el espacio de parámetros bivariado \((\theta,\mu)\) se puede factorizar como sigue:

\[p(\theta,\mu)=p(\theta|\mu)p(\mu)\]

Por lo tanto \[p(\theta,\mu|x)=\frac{p(x|\theta,\mu)p(\theta,\mu)}{p(x)}\] \[=\frac{p(x|\theta)p(\theta|\mu)p(\mu)}{p(x)}\] El siguiente modelo gráfico resume las independencias condicionales de la última ecuación:

library(DiagrammeR)
#> Warning: package 'DiagrammeR' was built under R version 3.5.1
grViz("
digraph betabinom {

  # a 'graph' statement
  graph [overlap = true]

  node [shape = box,penwidth = 0.1, fixedsize = true, fontsize = 9,
        fontname = Helvetica, width = 0.2, height = 0.2]
  AB[label = 'A,B']; 

  node [shape = circle,
        fixedsize = true, fontsize = 9,
        fontname = Helvetica, width = 0.2] // sets as circles
  mu[label = 'μ']; theta[label = 'θ']; x;

  # several 'edge' statements
  edge[color = grey, arrowsize = 0.5, penwidth = 0.5]
  AB->mu; mu->theta; theta->x;
}
")

Aproximación por cuadrícula

En el caso jerárquico, no se puede derivar la distribución posterior de manera analítica pero si los parámetros e hiperparámetros toman un número finito de
valores y no hay muchos de ellos, podemos aproximar la posterior usando aproximación por cuadrícula.

A continuación graficamos las distribuciones correspondientes al caso en que la distribución inicial del hiperparámetro \(\mu\) tiene la forma de una distribución \(Beta(2, 2)\), es decir, creemos que la media de la acuñadora \(\mu\) es 0.5, pero existe bastante incertidumbre acerca del valor central.

La distribución de \(\theta\), esto es, la distribución inicial que refleja la dependencia entre \(\theta\) y \(\mu\), se expresa por medio de otra distribución beta; en el ejemplo utilizaremos el valor fijo K=100:

\[\theta|\mu \sim Beta(\mu 100, (1-\mu)100)\]

Esta distribución inicial expresa un alto grado de certeza que una acuñadora con hiperparámetro \(\mu\) genera monedas con sesgo cercano a \(\mu\):

A_mu <- 2; B_mu <- 2; K <- 100

# Distribución inicial conjunta p(theta, mu)
p_conjunta <- function(mu, theta, A_mu, B_mu, K){
  # marginal p(mu)
  p_mu <- dbeta(mu, A_mu, B_mu)
  # condicional p(theta | mu)
  p_theta_mu <- dbeta(theta, mu * K, (1 - mu) * K)
  p_mu * p_theta_mu
}

grid <- expand.grid(theta = seq(0.01, 0.99, 0.005), 
  mu = seq(0.01, 0.99, 0.005))

grid_inicial <- grid %>% 
  mutate(p_inicial = p_conjunta(mu, theta, A_mu, B_mu, K))

plot_conj <- ggplot(grid_inicial, aes(x = theta, y = mu, z = p_inicial)) + 
  stat_contour(binwidth = 2, aes(color = ..level..)) +
  scale_x_continuous(expression(theta), limits = c(0, 1)) +
  scale_y_continuous(expression(mu), limits = c(0, 1)) +
  scale_color_gradient(expression(p(theta,mu)), guide = FALSE) 

grid_mu <- grid_inicial %>%
  group_by(mu) %>%
  summarise(p_mu = sum(p_inicial)) %>%
  ungroup() %>%
  mutate(p_mu = p_mu / sum(p_mu))

plot_mu <- ggplot(grid_mu, aes(x = mu, y = p_mu)) +
  geom_path() +
  labs(y = expression(p(mu)), 
       x = expression(mu)) + 
  ylim(0, 0.015) + coord_flip()

mus <- rbeta(5000, A_mu, B_mu)
sims_marg <- data.frame(sims = 1:5000, 
  sims_marg = rbeta(5000, mus * 100, (1 - mus) * 100))

plot_theta <- ggplot(sims_marg, aes(x = sims_marg)) + 
  geom_density(adjust = 2) +
  labs(y = expression(p(theta)), 
       x = expression(theta))

# p(theta|mu=0.75)
sims_cond_1 <- rbeta(5000, 0.75 * 100, (1 - 0.75) * 100)
# p(theta|mu=0.25)
sims_cond_2 <- rbeta(5000, 0.25 * 100, (1 - 0.25) * 100)

# marginal = sims_marg

sims <- data.frame(sim = 1:5000, 
  cond_1 = sims_cond_1, 
  cond_2 = sims_cond_2) %>%
  gather(dist, values, -sim) %>%
  mutate(mu = ifelse(dist == "cond_1", 0.75, 0.25))


plots_marginales <- ggplot(sims, aes(x = values)) +
  labs(y = expression(p(paste(theta, l, mu))), 
       x = expression(theta)) + 
  geom_density(adjust = 2) + 
  facet_wrap(~ mu, ncol = 1)

grid.arrange(plot_conj, plot_mu, plot_theta, plots_marginales, ncol=2)

Verosimilitud.


likeBern <- function(z, N){
  function(theta){
    theta ^ z * (1 - theta) ^ (N - z)
  }
}

# Valores observados
z <- 9; N <- 12
mi_like <- likeBern(z, N)

grid_like <- expand.grid(x = seq(0.01, 1, 0.05), y = seq(0.01, 1, 0.05))
grid_like <- grid_like %>% 
  mutate(z = mi_like(x))

ggplot(grid_like, aes(x = x, y = y, z = z)) + 
  stat_contour(aes(color = ..level..)) +
  scale_x_continuous(expression(theta), limits = c(0, 1)) +
  scale_y_continuous(expression(mu), limits = c(0, 1)) +
  scale_color_gradient(expression(L(theta,mu)), guide = FALSE) +
  coord_fixed()


# cuadrícula
grid_post <- grid_inicial %>%
  mutate(
    prior = p_inicial / sum(p_inicial),
    Like = theta ^ z * (1 - theta) ^ (N - z), # verosimilitud 
    posterior = (Like * prior) / sum(Like * prior)
  )  

post_conj <- ggplot(grid_post, aes(x = theta, y = mu, z = posterior)) + 
  stat_contour(aes(color = ..level..)) +
  scale_x_continuous(expression(theta), limits = c(0, 1)) +
  scale_y_continuous(expression(mu), limits = c(0, 1)) +
  scale_color_gradient(expression(p(theta,mu)), guide = FALSE)

Posterior.

# head(grid_post)

grid_mu <- grid_post %>%
  group_by(mu) %>%
  summarise(post_mu = sum(posterior)) %>%
  ungroup() %>%
  mutate(post_mu = post_mu / sum(post_mu))

post_mu <- ggplot(grid_mu, aes(x = mu, y = post_mu)) +
  geom_path() +
  labs(y = expression(p(paste(mu, l, x))), 
       x = expression(mu))+ coord_flip()

grid_theta <- grid_post %>%
  group_by(theta) %>%
  summarise(post_theta = sum(posterior)) %>%
  ungroup() %>%
  mutate(post_theta = post_theta / sum(post_theta))

post_theta <- ggplot(grid_theta, aes(x = theta, y = post_theta)) +
  geom_path() +
  labs(y = expression(p(paste(theta, l, x))), 
       x = expression(theta))

grid_theta_m <- grid_post %>%
  filter(mu == 0.75 | mu == 0.25) %>%
  group_by(theta, mu) %>%
  summarise(post_theta = sum(posterior)) %>%
  group_by(mu) %>%
  mutate(post_theta = post_theta / sum(post_theta))

post_marginales <- ggplot(grid_theta_m, aes(x = theta, y = post_theta)) +
  geom_path() + 
  facet_wrap(~ mu, ncol = 1) +
  labs(y = expression(p(paste(theta, l, x, mu))), 
       x = expression(theta))

Ejemplo: Múltiples monedas de una misma fábrica

La sección anterior considera el escenario en que lanzamos una moneda y hacemos inferencia del parámetro de sesgo \(\theta\) y del hiperparámetro \(\mu\). Ahora consideramos recolectar información de múltiples monedas, si cada moneda tiene su propio sesgo \(\theta_j\) entonces debemos estimar un parámetro distinto para cada moneda.

Suponemos que todas las monedas provienen de la misma fábrica, esto implica que tenemos la mismas creencias iniciales \(\mu\) para todas las monedas. Suponemos también que cada moneda se acuña de manera independiente, esto es, que condicional al parámetro \(\mu\) los parámetros \(\theta_j\) son independientes en nuestras creencias iniciales.

Aproximación por cuadrícula

Supongamos que tenemos dos monedas de la misma fábrica. El objetivo es estimar los sesgos \(\theta_1\), \(\theta_2\) de las dos monedas y estimar simultáneamente el parámetro \(\mu\) correspondiente a la casa de moneda que las fabricó.

Inicial. Asumimos la misma distribución inicial para \(\theta_1\) y \(\theta_2\).

Verosimilitud.


z_1 <- 3; N_1 <- 15
z_2 <- 4; N_2 <- 5

mi_like_1 <- likeBern(z_1, N_1)
mi_like_2 <- likeBern(z_2, N_2)

grid_like <- grid_like %>% 
  mutate(L_1 = mi_like_1(x), 
         L_2 = mi_like_2(x))

plot_like_1 <- ggplot(grid_like, aes(x = x, y = y, z = L_1)) + 
  stat_contour(aes(color = ..level..)) +
  scale_x_continuous(expression(theta[1]), limits = c(0, 1)) +
  scale_y_continuous(expression(mu), limits = c(0, 1)) +
  scale_color_gradient(expression(L(theta,mu)), guide = FALSE)

plot_like_2 <- ggplot(grid_like, aes(x = x, y = y, z = L_2)) + 
  stat_contour(aes(color = ..level..)) +
  scale_x_continuous(expression(theta[1]), limits = c(0, 1)) +
  scale_y_continuous(expression(mu), limits = c(0, 1)) +
  scale_color_gradient(expression(L(theta,mu)), guide = FALSE)

grid.arrange(plot_like_1, plot_like_2, ncol = 3)

Posterior.

JAGS

La sección anterior utilizó un modelo simplificado con el objetivo de poder visualizar el espacio de parámetros. Ahora incluiremos más parámetros con el objetivo de hacer el problema más realista. En los ejemplos anteriores fijamos el grado de dependencia entre \(\mu\) y \(\theta\) de manera arbitraria, a través de \(K\), de tal manera que si \(K\) era grande, los valores \(\theta_j\) individuales se situaban cerca de \(\mu\) y cuando \(K\) era pequeña se permitía más variación.

En situaciones reales no conocemos el valor de \(K\) por adelantado, por lo que dejamos que los datos nos informen acerca de valores creíbles para \(K\). Intuitivamente, cuando la proporción de águilas observadas es similar a lo largo de las monedas, tenemos evidencia de que \(K\) es grande, mientras que si estas proporciones difieren mucho, tenemos evidencia de que \(K\) es pequeña. Debido a que \(K\) pasará de ser una constante a ser un parámetro lo llamaremos \(\kappa\).

La distribución inicial para \(\kappa\) será una Gamma. Elegimos \(\kappa \sim Gamma(1, 0.1)\); ésta tiene media 10 y desviación estándar 10.

library(R2jags)
#> Warning: package 'R2jags' was built under R version 3.5.1
#> Loading required package: rjags
#> Loading required package: coda
#> Linked to JAGS 4.3.0
#> Loaded modules: basemod,bugs
#> 
#> Attaching package: 'R2jags'
#> The following object is masked from 'package:coda':
#> 
#>     traceplot

modelo_jer.txt <-
'
model{
  for(t in 1:N){
    x[t] ~ dbern(theta[coin[t]])
  }
  for(j in 1:nCoins){
    theta[j] ~ dbeta(a, b)
  }
  a <- mu * kappa
  b <- (1 - mu) * kappa
  # A_mu = 2, B_mu = 2
  mu ~ dbeta(2, 2)
  kappa ~ dgamma(1, 0.1)
}
'

Los datos consisten en 5 monedas, cada una se lanza 5 veces cada una, resultando 4 de ellas en 1 águila y 4 soles y otra en 3 águilas y 2 soles.

cat(modelo_jer.txt, file = 'modelo_jer.bugs')

x <- c(0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 1, 1)
coin <- c(rep(1, 5), rep(2, 5), rep(3, 5), rep(4, 5), rep(5, 5))

jags.inits <- function(){
  list("mu" = runif(1, 0.1, 0.9),
    "kappa" = runif(1, 5, 20))
}

jags_fit <- jags(
  model.file = "modelo_jer.bugs",    # modelo de JAGS
  inits = jags.inits,   # valores iniciales
  data = list(x = x, coin = coin, nCoins = 5,  N = 25),    # lista con los datos
  parameters.to.save = c("mu", "kappa", "theta"),  # parámetros por guardar
  n.chains = 3,   # número de cadenas
  n.iter = 10000,    # número de pasos
  n.burnin = 1000   # calentamiento de la cadena
  )
#> module glm loaded
#> Compiling model graph
#>    Resolving undeclared variables
#>    Allocating nodes
#> Graph information:
#>    Observed stochastic nodes: 25
#>    Unobserved stochastic nodes: 7
#>    Total graph size: 65
#> 
#> Initializing model

traceplot(jags_fit, varname = c("kappa", "mu", "theta"))


jags_fit
#> Inference for Bugs model at "modelo_jer.bugs", fit using jags,
#>  3 chains, each with 10000 iterations (first 1000 discarded), n.thin = 9
#>  n.sims = 3000 iterations saved
#>          mu.vect sd.vect   2.5%    25%    50%    75%  97.5% Rhat n.eff
#> kappa     14.874  10.838  2.213  6.824 12.164 19.819 42.307    1  1500
#> mu         0.335   0.098  0.163  0.267  0.328  0.400  0.537    1  3000
#> theta[1]   0.287   0.128  0.070  0.195  0.277  0.369  0.560    1  3000
#> theta[2]   0.291   0.129  0.067  0.195  0.283  0.377  0.559    1  3000
#> theta[3]   0.290   0.130  0.069  0.196  0.281  0.373  0.573    1  3000
#> theta[4]   0.292   0.127  0.073  0.200  0.281  0.377  0.562    1  3000
#> theta[5]   0.419   0.144  0.165  0.319  0.408  0.512  0.733    1  3000
#> deviance  30.354   2.124 27.481 28.873 29.940 31.399 35.591    1  3000
#> 
#> For each parameter, n.eff is a crude measure of effective sample size,
#> and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
#> 
#> DIC info (using the rule, pD = var(deviance)/2)
#> pD = 2.3 and DIC = 32.6
#> DIC is an estimate of expected predictive error (lower deviance is better).

Podemos hacer histogramas de las distribuciones posteriores de los parámetros:

sims_df <- data.frame(n_sim = 1:jags_fit$BUGSoutput$n.sims,
  jags_fit$BUGSoutput$sims.matrix) %>% 
  dplyr::select(-deviance) %>%
  gather(parametro, value, -n_sim)

ggplot(filter(sims_df, parametro != "kappa"), aes(x = value)) +
  geom_histogram(alpha = 0.8) + facet_wrap(~ parametro)
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.



ggplot(filter(sims_df, parametro == "kappa"), aes(x = value)) +
  geom_histogram(alpha = 0.8)
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Ejemplo: Función de Verosimilitud Poisson

options(digits = 3)
library(ggplot2)
library(plyr)
library(dplyr)
library(reshape2)
library(tidyr)
library(gridExtra)
library(MASS)

En esta sección veremos un problema de estimación de tasas de mortalidad.

Plantearemos 3 alternativas de modelación para resolver el problema: modelo de unidades iguales, modelo de unidades independientes y, finalmente, modelo jerárquico.

La base de datos tiene información de todas las cirugías de transplante de corazón llevadas a cabo en Estados Unidos en un periodo de 24 meses, entre octubre de 1987 y diciembre de 1989.

La forma en como se construyó la base de datos es como sigue:

library(LearnBayes)
data(hearttransplants)
heart <- hearttransplants
head(heart)
#>      e y
#> 1  532 0
#> 2  584 0
#> 3  672 2
#> 4  722 1
#> 5  904 1
#> 6 1236 0

El objetivo es obtener buenas estimaciones de las tasas de mortalidad de cada hospital. Antes de comenzar a ajustar modelos complejos vale la pena observar los datos. El cociente \(\{y_{j}/e_{j}\}\) es el número observado de muertes por unidad de exposición y se puede ver como una estimación puntual de la tasa de mortalidad.

La siguiente figura grafica, en el eje vertical el cociente \(\{y_{j}/e_{j}\}\) multiplicado por 1000 (con la intención de que la tasa indique número de muertes por 1000 expuestos), y en el eje horizontal el número de expuestos \(\{e_{j}\}\) -en escala logarítmica- para los 94 hospitales. Cada punto representa un hospital y esta etiquetado con el número de muertes observadas \(\{y_{j}\}\).

ggplot(heart, aes(x = log(e), y = 1000* y / e, color = factor(y), label = y)) +
  geom_text(show.legend = FALSE)  +
  scale_x_continuous("Número de expuestos (e)", labels = exp, 
    breaks = c(log(700), log(700*2), log(700*2^2), log(700*2^3), 
      log(700*2^4)))

En la gráfica podemos notar que las tasas estimadas son muy variables, especialmente para hospitales con baja exposición. Observemos también que la mayoría de los hospitales que no experimentaron muertes tienen bajo número de expuestos.

La variabilidad de las tasas de mortalidad y los hospitales sin muertes sugieren un problema de tamaño de muestra: Consideremos un hospital con 700 expuestos. La muerte por transplante de corazón no es común, el promedio nacional es de 9 por cada 10,000 expuestos, por lo que con 700 expuestos es probable que no se presente ninguna muerte, en este caso el hospital pertenece al 10% de los hospitales con menor tasa de mortalidad. Sin embargo, existe la posibilidad de que se observe una muerte, con lo cual el hospital tendrá una tasa de mortalidad que es lo suficientemente alta para que quede en el 25% de los hospitales con mayor tasa de mortalidad observada.

Una vez reconocido el problema de utilizar los datos crudos para estimar las tasas de mortalidad plantearemos 3 alternativas de modelación bayesiana. En todos los casos modelaremos cada observación del número de muertes \(y_i\) mediante una distribución \(f(y_i|\lambda_i)\) con parámetro \(\lambda_i\) que representa la tasa de mortalidad. El objetivo es estimar estas tasas de mortalidad.

  1. Modelo de unidades iguales. Considera que los estudios son replicas idénticas una de otra, en este sentido vemos a las observaciones de todos los estudios como muestras independientes de una misma población y consecuentemente todas las observaciones tendrán la misma tasa de mortalidad \(\lambda_j = \lambda\).

\[ \begin{eqnarray} \nonumber y_{j} &\sim& f(y_{j}|\lambda),\\ \nonumber \lambda &\sim& g(\lambda). \end{eqnarray} \] donde \(g(\lambda)\) es la distribución inicial del único parámetro \(\lambda\).

El modelo gráfico correspondiente a este enfoque (omitiendo la distribución de \(\lambda\)) es:

  1. Modelo de unidades independientes. El extremo opuesto al modelo de unidades iguales considera que los estudios son tan diferentes que los resultados de cada uno no proveen información acerca de los resultados de ningún otro y, por tanto, realizamos estimaciones individuales para cada \(\lambda_{j}\), \[ \begin{align} \nonumber y_{j} &\sim f(y_{j}|\lambda_{j}),\\ \nonumber \lambda_j &\stackrel{\mathrm{iid}}{\sim} g(\lambda_j). \end{align} \] donde \(g(\lambda_j)\) es la distribución inicial del parámetro \(\lambda_j\). Note que asumimos la misma distribución inicial para todos los parámetros.

El modelo gráfico correspondiente (omitiendo las distribución de las \(\lambda_j\)) es:

  1. Modelo jerárquico. Es un análisis intermedio que considera los parámetros de interés \(\lambda_{j}\) como provenientes de una distribución común, \[ \begin{eqnarray} \nonumber y_{j} &\sim& f(y_{j}|\lambda_{j}),\\ \nonumber \lambda_{j} &\sim& g(\lambda_j|\theta),\\ \theta &\sim& h(\theta). \end{eqnarray} \] donde \(g(\lambda_j|\theta)\) es la distribución inicial del parámetro \(\lambda_j\) con parámetro \(\theta\) y \(h(\theta)\) es la distribución inicial de \(\theta\).

El modelo gráfico asociado (omitiendo la distribución de \(\theta\)) es:

De esta manera, la probabilidad conjunta del modelo refleja una dependencia entre los parámetros al mismo tiempo que permite variaciones entre los estudios. Como resultado se estima una \(\lambda_{j}\) diferente para cada estudio usando información de todos.

Es importante notar que la estructura del modelo jerárquico facilita la implementación del muestreador de Gibbs debido a que es inmediato factorizar la distribución conjunta en distribuciones marginales condicionales completas.

Modelo de unidades iguales

Supongamos que las tasas de mortalidad son iguales a lo largo de los hospitales. Estimaremos la tasa de mortalidad con el siguiente modelo para los datos, \[ \begin{eqnarray} \nonumber y_{j}|\lambda \sim Poisson(e_{j}\lambda), \end{eqnarray} \]

donde \(y_{j}\) es el número de muertes observadas por transplante corazón en el hospital \(j\), \(e_{j}\) es el número de expuestos, y \(\lambda\) es la tasa de mortalidad, medida en número de muertes por unidad de exposción y común para todos los hospitales. La distribución Poisson tiene media \(e_j \lambda\).

Sea \(y=(y_1,...,y_{94})\), la función de verosimilitud es entonces,

\[ \begin{eqnarray} \nonumber \mathcal{L}(\lambda) \equiv f(y|\lambda) = \prod_{i=1}^{94}f(y_{i}|\lambda) = \prod_{i=1}^{94}\bigg(\frac{exp(-e_{i}\lambda)(e_{i}\lambda)^{y_{i}}}{y_{i}!}\bigg) \end{eqnarray} \]

Debido a que no contamos con información inicial acerca de la tasa de mortalidad, asignamos a \(\lambda\) una distribución inicial no informativa, de la forma \[g(\lambda)\propto\frac{1}{\lambda}.\]

Usamos el Teorema de Bayes para calcular la densidad posterior de \(\lambda\),

\[ \begin{eqnarray} \nonumber g(\lambda|y) &\propto& g(\lambda)f(y|\lambda) \\ \nonumber &=& \frac{1}{\lambda}\prod_{i=1}^{94}\bigg(\frac{exp(-e_{i}\lambda)(e_{i}\lambda)^{y_{i}}}{y_{i}!}\bigg)\\ \nonumber &\propto& \lambda^{(\sum_{i=1}^{94}y_{j}-1)}exp\bigg(-\sum_{i=1}^{94}e_{j}\lambda\bigg), \end{eqnarray} \]

Con el núcleo de esta función, identificamos la distribución posterior como una \(Gamma(\sum_{i=1}^{94}y_{i},\sum_{i=1}^{94}e_{i})\), donde expresamos la función de densidad de una distribución \(Gamma(\alpha, \beta)\) usando el parámetro de forma \(\alpha\) y el inverso del parámetro de escala \(\beta\), de manera que la función de densidad de probabilidad posterior es,

\[ \begin{align} \nonumber &f(x|\alpha, \beta) = \beta^{\alpha}\frac{1}{\Gamma(\alpha)}x^{\alpha-1}e^{-\beta x}I_{(0,\infty)}(x)\\ \nonumber &\mbox{para }x \geq 0 \mbox{ y }\alpha \mbox{, }\beta > 0. \end{align} \] donde \(I(x)\) es la función indicadora tal que \[ I_{(0,\infty)}(x) = \left\{ \begin{align} \nonumber & 1 \quad \text{si} \quad 0<x<\infty \\ \nonumber & 0 \quad e.o.c. \end{align} \right. \]

Predicción:

Para verificar el ajuste del modelo utilizaremos la distribución predictiva posterior. Si \(y_{j}^*\) denota el número de muertes por transplante de corazón en el hospital \(j\) con exposición \(e_{j}\) en una muestra futura, la distribución predictiva posterior es

\[ \begin{eqnarray} \nonumber f(y_{j}^*|y) = \int f(y_{j}^*|\lambda)g(\lambda|y)d\lambda, \end{eqnarray} \]

donde \(f(y_j^*|\lambda) \sim Poisson(e_i \lambda)\).

La distribución predictiva posterior \(f(y_{j}^*|y)\) se puede interpretar como la función de verosimilitud de observaciones futuras basadas en el modelo ajustado (en el lugar de los parámetros del modelo). En este ejemplo, la densidad \(f(y_{j}^*|y)\) representa el número de transplantes que se predecirían para un hospital con exposición \(e_{j}\). Entonces, si el número observado de muertes \(y_{j}\) no está en las colas de la distribución predictiva posterior, diríamos que la observación es consistente con el modelo observado (es lo que predeciría el modelo ajustado).

# La densidad posterior es Gamma con los siguientes parámetros
c(sum(heart$y), sum(heart$e))
#> [1]    277 294681
heart$hospital <- 1:94
# Simulamos 1000 muestras de la densidad posterior de lambda
lambdas <- rgamma(1000, shape = sum(heart$y), rate = sum(heart$e))
# ahora para cada hospital simulamos muestras de una distribución Poisson
# con media e_i*lambda
sims <- ddply(heart, "hospital", transform, sims = rpois(1000, e*lambdas))
# tomamos una muestra de 10 hospitales
set.seed(242)
hosps <- sample(1:94, 10)
sims.2 <- subset(sims, hospital %in% hosps)
heart.2 <- subset(heart, hospital %in% hosps) 
head(sims.2)
#>     e y hospital sims
#> 1 532 0        1    1
#> 2 532 0        1    0
#> 3 532 0        1    0
#> 4 532 0        1    1
#> 5 532 0        1    0
#> 6 532 0        1    0

ggplot(sims.2, aes(x = sims)) + 
  geom_histogram(aes(y = ..density..), binwidth = 1, color = "darkgray", 
    fill = "darkgray") +
  facet_wrap(~ hospital, nrow = 2) +
  geom_segment(data = heart.2, aes(x = y, xend = y, y = 0, yend = 0.5), 
    color = "red") + 
  geom_text(data = heart.2, aes(x = 10, y = 0.4, label = paste("e =", e)), 
    size = 2.7)

La figura muestra los histogramas, obtenidos con simulación, de la distribución predictiva posterior de una muestra de 10 hospitales de los 94. Para cada hospital, se escribió el número de expuestos, \(e\), y se grafica una línea vertical en el número de muertes observado. Notemos que en muchos de los histogramas el número de muertes observado se encuentra en las colas de las distribuciones, lo que sugiere que nuestras observaciones son inconsistentes con el modelo ajustado.

Ahora examinaremos la consistencia de las muertes observadas, \(y_{j}\), para todos los hospitales. Para cada distribución predictiva posterior calculamos la probabilidad de que una observación futura \(y_{j}^*\) sea al menos tan extrema como \(y_{j}\), estas probabilidades son comúnmente llamadas valores-p predictivos:

\[ \begin{eqnarray} \nonumber P(extremos) = min\{P(y_{j}^* \leq y_{j}),P(y_{j}^*\geq y_{j})\} \end{eqnarray} \]

y graficamos estos valores-p para cada uno de los 94 hospitales:

# Para calcular las p predictiva podemos usar las simulaciones
p.pred <- ddply(sims, "hospital", summarise, p = min(sum(sims <= y) / 1000, 
  sum(sims >= y) / 1000))
head(p.pred)
#>   hospital     p
#> 1        1 0.598
#> 2        2 0.586
#> 3        3 0.146
#> 4        4 0.499
#> 5        5 0.560
#> 6        6 0.334
p.pred.2 <- join(heart, p.pred, by = "hospital")
ggplot(p.pred.2, aes(x = log(e), y = p)) + 
  geom_point() +
  scale_x_continuous("Número de expuestos (e)", labels = exp, 
    breaks = c(log(700), log(700*2), log(700*2^2), log(700*2^3), 
      log(700*2^4))) + 
  ylab("P(extremos)") +
  geom_hline(yintercept = .15, colour = "red", size = 0.4) +
  ylim(0, .7)

En la figura anterior graficamos las probabilidades de extremos (calculadas con simulación) en el eje vertical, contra el número de exposición en escala logarítmica de cada hospital en el eje horizontal. Notemos que muchas de estas probabilidades son pequeñas, el 28% son menores a 0.15, lo que indica que para el 28% de los hospitales el número de muertes observado \(y_{j}\) está en la cola de la distribución y por consiguiente el modelo es inadecuado.

Modelo de unidades independientes

Consideremos ahora que los hospitales son independientes; estimaremos la tasa de mortalidad para cada hospital con el modelo,

\[ \begin{eqnarray} \nonumber y_{j}|\lambda_j \sim Poisson(e_{j}\lambda_{j}), \end{eqnarray} \]

donde \(y_{j}\) es el número de muertes observadas por transplante corazón en el hospital \(j\), \(e_{j}\) es el número de expuestos, y \(\lambda_{j}\) es la tasa de mortalidad, medida en número de muertes por unidad de exposición, con \(j=1,...,94\).

Note que, ya que se trata de unidades independientes, ajustaremos un modelo para cada subíndice \(j\). Por lo tanto, la función de verosimilitud para el parámetro \(\lambda_j\) coincide con la distribución de la observación \(y_j\): \[ \mathcal{L}(\lambda_j) = f(y_j|\lambda_j) = Poisson(e_j \lambda_j)\]

Por el momento asignamos a \(\lambda_{j}\) una distribución inicial \(Gamma(\alpha_{0}, \beta_{0})\) que es conjugada de la distribución Poisson,

\[ \begin{eqnarray} \nonumber g(\lambda_{j}|\alpha_{0},\beta_{0}) = \frac{\beta_{0}^{\alpha_{0}} \lambda_{j}^{\alpha_{0}-1}exp(-\lambda_{j} \beta_{0})}{\Gamma(\alpha_{0})}, \quad \lambda_{j}>0. \end{eqnarray} \]

Calculamos la distribución posterior de \(\lambda_{j}\) usando el Teorema de Bayes,

\[ \begin{eqnarray} \nonumber g(\lambda_{j}|y_{j},\alpha_{0},\beta_{0}) &\propto& g(\lambda_{j}|\alpha_{0},\beta_{0})f(y_{j}|\lambda_{j})\\ \nonumber &=& \frac{\beta_{0}^{\alpha_{0}} \lambda_{j}^{\alpha_{0}-1}exp(-\lambda_{j} \beta_{0})}{\Gamma(\alpha_{0})}\frac{exp(-e_{j}\lambda_{j})(e_{j}\lambda_{j})^{y_{j}}}{y_{j}!}\\ \nonumber &\propto& \lambda_{j}^{\alpha_{0}+y_{j}-1}exp(-\lambda_{j}(\beta_{0}+e_{j})), \end{eqnarray} \]

y obtenemos que la distribución posterior es \(\lambda_{j}|y_{j} \sim Gamma(\alpha_{0}+y_{j}, \beta_{0}+e_{j})\).

La media de la distribución posterior es

\[ \begin{eqnarray} \nonumber E(\lambda_{j}|y_{j},\alpha_{0},\beta_{0}) &=& \frac{\alpha_{0}+y_{j}}{\beta_{0}+e_{j}}\\ \label{eqn:pond.indep} &=& (1-A_{j})\frac{y_{j}}{e_{j}}+A_{j}\frac{\alpha_{0}}{\beta_{0}} \end{eqnarray} \]

donde

\[ \begin{eqnarray} A_{j}=\frac{\beta_{0}}{\beta_{0}+e_{j}}, \end{eqnarray} \] y su varianza es

\[ \begin{eqnarray} \nonumber Var(\lambda_{j} |y_{j},\alpha_{0},\beta_{0}) = \frac{\alpha_{0}+y_{j}}{(\beta_{0}+e_{j})^2}. \end{eqnarray} \]

Notemos que podemos escribir la media posterior como un promedio ponderado de la tasa observada \(y_{j}/e_{j}\) y la media de la distribución inicial \(\alpha_{0} / \beta_{0}\). El factor \(A_{j}\) se puede interpretar como el encogimiento hacia la media inicial y depende del número de exposición \(e_j\) y del parámetro de escala \(\beta_{0}\).

Efecto de \({\beta_{0}}\) y \({e_{j}}\) en la distribución posterior

Consideremos la media y varianza posteriores de \(\lambda_{j}\) para un hospital particular. Tomando \(\alpha_{0}\) fija, valores mayores de \(\beta_{0}\) corresponden a una menor varianza en la distribución inicial ya que \[Var(\lambda_{j} | \alpha_{0},\beta_{0}) = \alpha_{0}/\beta_{0}^2\]

La varianza de la distribución inicial se refleja en la media de la distribución posterior de \(\lambda_{j}\) a través de \(\beta_{0}\): menor varianza (mayor \(\beta_{0}\)) corresponde a mayor encogimineto hacia la media inicial. Este efecto de \(\beta_{0}\) concuerda con la incertidumbre de nuestro conocimiento, pues menor varianza en la distribución inicial indica menor incertidumbre y la aportación de la media inicial a la media posterior es más importante.

La elección de \(\beta_{0}\) también afecta la varianza, pues un menor valor de \(\beta_{0}\) implica una mayor varianza tanto en la distribución inicial como en la distribución posterior.

Por su parte, tomando \(\alpha_{0}\) y \(\beta_{0}\) fijas, el efecto del número de expuestos \(e_{j}\) sobre la media de la distribución posterior de \(\lambda_{j}\) actúa de manera contraria a \(\beta_{0}\). Un mayor número de expuestos tiene como consecuencia un menor encogimiento hacia la media inicial, dando mayor importancia a la tasa observada \(y_{j}/e_{j}\). Esto refleja que más expuestos implican más información proveniente de la muestra, restando importancia a nuestro conocimiento inicial.

En cuanto a la varianza de la distribución posterior, es inversamente proporcional al número de expuestos, indicando nuevamente que más expuestos implica más conocimiento de la muestra y por tanto menor incertidumbre.

En la siguiente fugura mostraremos el encogimiento hacia la media de la distribución inicial bajo distintos escenarios de \(\beta_{0}\), cada punto representa un hospital y el color indica a qué valor de \(\beta_{0}\) corresponde.

encoge <- ddply(heart, "hospital", transform, enc.1 = 1000 / (1000 + e), 
  enc.2 = 100 / (100 + e), enc.3 = 10 / (10 + e), enc.4 = 0.01 / (0.01 + e))
head(encoge)
#>      e y hospital enc.1  enc.2   enc.3    enc.4
#> 1  532 0        1 0.653 0.1582 0.01845 1.88e-05
#> 2  584 0        2 0.631 0.1462 0.01684 1.71e-05
#> 3  672 2        3 0.598 0.1295 0.01466 1.49e-05
#> 4  722 1        4 0.581 0.1217 0.01366 1.39e-05
#> 5  904 1        5 0.525 0.0996 0.01094 1.11e-05
#> 6 1236 0        6 0.447 0.0749 0.00803 8.09e-06
encoge.m <- melt(encoge, id = c("e", "y", "hospital"))
head(encoge.m)
#>      e y hospital variable value
#> 1  532 0        1    enc.1 0.653
#> 2  584 0        2    enc.1 0.631
#> 3  672 2        3    enc.1 0.598
#> 4  722 1        4    enc.1 0.581
#> 5  904 1        5    enc.1 0.525
#> 6 1236 0        6    enc.1 0.447

ggplot(encoge.m, aes(x = log(e), y = value, color = variable)) +
  geom_point(alpha = 0.9) +
  scale_x_continuous("Número de expuestos (e)", labels = exp, 
    breaks = c(log(700), log(700*2), log(700*2^2), log(700*2^3), 
      log(700*2^4))) +
  scale_color_hue(expression(beta[0]), labels = c("1000", "100", "10", "0.01")) +
  ylab("encogimiento (A)")

En esta gráfica podemos apreciar el mayor encogimiento para valores mayores de \(\beta_{0}\) y el decaimiento en el encogimiento conforme aumenta el número de expuestos.

Distribución inicial

Establecimos que por facilidad se utilizará una inicial \(Gamma(\alpha_{0},\beta_{0})\), sin embargo hace falta asignar valores a los parámetros iniciales. Analizaremos 3 distintas combinaciones de parámetros iniciales y su efecto en las estimaciones posteriores.

Para cada combinación de parámetros la siguiente tabla muestra los deciles de 2000 simulaciones de una distribución \(Gamma(\alpha_0,\beta_{0})\) y los deciles de las muertes observadas por cada 1000 expuestos considerando todos los hospitales; multiplicamos las simulaciones por 1000 para que indiquen tasas de mortalidad medidas en número de muertes por 1000 expuestos.

El propósito de la tabla es describir la forma de la distribución \(Gamma\) al cambiar los parámetros, por ejemplo, una \(Gamma(0.01,0.01)\) es muy plana y podríamos describirla como poco informativa, mientras que una \(Gamma(1,1000)\) tiene una forma más cercana a las tasas de mortalidad observadas.

decil \((0.01,0.01)\) \((0.5,0.01)\) \((1,1000)\) Observados
min 0 0 0 0
10 0 890.3 0.1 0
20 0 3134.7 0.2 0.3
30 0 7769.0 0.3 0.5
40 0 15036.0 0.5 0.7
50 0 24293.9 0.7 0.8
60 0 37306.0 0.9 1.1
70 0 57640.2 1.2 1.4
80 0 84255.5 1.6 1.7
90 4.2 140721.9 2.3 2.0
max 287307.6 592073.2 7.0 3.9

Ahora realizamos una gráfica para cada una de las 3 parejas de parámetros
\((\alpha_0,\beta_{0})\) donde mostramos en negro las tasas de mortalidad observadas, y en color las estimaciones posteriores de las tasas de mortalidad con intervalos del 95% de probabilidad, ambas medidas en número de muertes por cada 1000 expuestos.

Cada punto representa un hospital y el color corresponde al número de muertes observadas, \(\{y_{j}\}\). Regresaremos a esta gráfica más adelante pero por lo pronto observemos que tanto las estimaciones como los intervalos de probabilidad son muy diferentes al cambiar los parámetros de la distribución inicial.

# Procedemos como antes, para cada combinación de alfa y beta simulamos 1000 
# lambdas de la posterior
lambdas <- ddply(heart, "hospital", transform, 
  sims =  rgamma(n = 2000, shape = 0.01 + y, rate = 0.01 + e))
# Creamos intervalos con las simulaciones
int.post <- ddply(lambdas, "hospital", summarise, 
  int.izq = quantile(1000*sims, 0.025), int.der = quantile(1000*sims, 0.975), 
  media = 1000*(0.01 + y[1])/(0.01 + e[1]))
int.post$comb <- "alpha = 0.01 beta = 0.01"
heart.int <- join(heart, int.post, by = "hospital")
ggplot(heart.int, aes(x = log(e), y = media, color = factor(y))) +
  geom_point() +
  geom_segment(aes(x = log(e), xend = log(e), y = int.izq, yend = int.der)) +
  geom_point(data = heart, aes(x = log(e), y = 1000*y/e), color = "black", 
    alpha = 0.6)


# alpha = 0.5, beta = 0.01
lambdas <- ddply(heart, "hospital", transform, 
  sims =  rgamma(n = 2000, shape = 0.5 + y, rate = 0.01 + e))
# Creamos intervalos con las simulaciones
int.post.2 <- ddply(lambdas, "hospital", summarise, 
  int.izq = quantile(1000*sims, 0.025), int.der = quantile(1000*sims, 0.975), 
  media = 1000*(0.5 + y[1])/(0.01 + e[1]))
int.post.2$comb <- "alpha = 0.5 beta = 0.01"

# alpha = 1, beta = 1000
lambdas <- ddply(heart, "hospital", transform, 
  sims =  rgamma(n = 2000, shape = 1 + y, rate = 1000 + e))
# Creamos intervalos con las simulaciones
int.post.3 <- ddply(lambdas, "hospital", summarise, 
  int.izq = quantile(1000*sims, 0.025), int.der = quantile(1000*sims, 0.975), 
  media = 1000*(0.01 + y[1])/(0.01 + e[1]))
int.post.3$comb <- "alpha = 1 beta = 1000"

int.post <- rbind(int.post, int.post.2, int.post.3)
heart.int <- join(heart, int.post, by = "hospital")
head(heart.int)
#>     e y hospital   int.izq int.der  media                     comb
#> 1 532 0        1 1.06e-145   0.113 0.0188 alpha = 0.01 beta = 0.01
#> 2 532 0        1  6.35e-04   4.474 0.9398  alpha = 0.5 beta = 0.01
#> 3 532 0        1  1.27e-02   2.386 0.0188    alpha = 1 beta = 1000
#> 4 584 0        2 1.39e-149   0.120 0.0171 alpha = 0.01 beta = 0.01
#> 5 584 0        2  1.03e-03   4.264 0.8561  alpha = 0.5 beta = 0.01
#> 6 584 0        2  1.84e-02   2.294 0.0171    alpha = 1 beta = 1000
ggplot(heart.int, aes(x = log(e), y = media, color = factor(y))) +
  geom_point() +
  geom_segment(aes(x = log(e), xend = log(e), y = int.izq, yend = int.der)) +
  facet_wrap(~ comb, nrow = 3) +
  geom_point(data = heart, aes(x = log(e), y = 1000*y/e), color = "black", 
    alpha = 0.6) +
  scale_x_continuous("Número de expuestos (e)", labels = exp, 
    breaks = c(log(700), log(700*2), log(700*2^2), log(700*2^3), 
      log(700*2^4))) + 
  scale_colour_hue("muertes obs. (y)") + 
  ylab("Muertes por 1000 expuestos")

Ahora analizaremos la elección de las parejas de parámetros iniciales. La única información con la que contamos para definir una distribución inicial es que la tasa de mortalidad por transplante de corazón debe ser positiva y no muy grande.

No informativas

Asignamos a \(\lambda_{j}\) una distribución inicial \(Gamma(0.01,0.01)\). La tabla de deciles indica que es una distribución muy plana y si observamos la gráfica de encogimiento notamos que para una inicial con este valor en el parámetro de escala \(\beta_{0}\), el encogimiento de la media posterior hacia la media inicial es muy cercano a cero, dando poca importancia a la media inicial \(\alpha_{0}/\beta_{0}=1\). Como consecuencia, las estimaciones posteriores de \(\lambda_{j}\) son casi iguales a las tasas observadas \(\{y_{j}/e_{j}\}\), y se presentan los problemas de tamaño de muestra notados al usar las tasas observadas como estimaciones de las tasas de mortalidad. Además, los intervalos de confianza para los hospitales que no experimentaron muertes son poco creíbles pues son mucho más chicos que el resto.

Semi-informativas

La siguiente distribución inicial que consideramos es una \(Gamma(0.5,0.01)\). En este caso, la elección de los parámetros se basó en la distribución inicial no informativa de Jeffreys, consiste en una \(Gamma(0.5,0)\). Es una distribución impropia por lo que cambiamos el parámetro de escala por \(0.01\) para obtener una inicial propia con varianza grande. Los resultados no son muy razonables, pues la media inicial es 50, mucho mayor a las tasas observada para todos los hospitales, provocando que las estimaciones del modelo sean mayores a las tasas observadas en todos los hospitales. Este efecto es contrario al que buscábamos al usar una distribución inicial no informativa pues la distribución inicial tiene un impacto muy fuerte en las estimaciones posteriores.

Informativas

Finalmente, consideramos una distribución inicial \(Gamma(1,1000)\), ésta inicial es informativa. Para obtener sus parámetros igualamos la media y la varianza teóricas de la distribución \(Gamma\) con la media y la varianza observadas en el conjunto de tasas de mortalidad tomando en cuenta todos los hospitales. Las estimaciones de las tasas de mortalidad que obtenemos parecen razonables; sin embargo, especificar la distribución inicial con la muestra tiene problemas lógicos y prácticos:

  1. los datos se están usando 2 veces, primero la información de todos los hospitales se usa para especificar la distribución inicial, y después la información de cada hospital se usa para estimar su tasa de mortalidad \(\lambda_{j}\), lo que ocasiona sobreestimación.

  2. De acuerdo a la lógica bayesiana no tiene sentido estimar \(\alpha_{0}\) y \(\beta_{0}\), pues son parte de las creencias iniciales y no deben depender de los datos.

Conclusión

Podemos concluir que el modelo de unidades independientes no es robusto para nuestros datos pues las estimaciones posteriores de las tasas de mortalidad son muy sensibles a la elección de los parámetros de la distribución inicial.

Por tanto, parece ser conveniente intentar mejorar las estimaciones individuales de \(\lambda_{j}\) usando la información de todos los hospitales. La manera correcta de hacerlo es establecer un modelo de probabilidad para todo el conjunto de parámetros
\(\{\alpha,\beta,\lambda_{1},...,\lambda_{94}\}\) y después realizar un análisis de la distribución conjunta.

Se llevará a cabo un análisis completo en la siguiente sección en donde se usará un modelo jerárquico.